/*-------------------------------------------------------------------------*/
/**
  @file		initial.c
  @author	Renaud Wanchoor
  @date	July 2004
  @version	$Revision: 2.9 $
  @brief	initial partitioning 
  
*/
/*--------------------------------------------------------------------------*/
#include "pagrid.h"
#include "execrefine_params.h"

/**
 * \brief Initial partitioning
 *
 * Find the processor assignment that minimizes the total execution time
 * of the application, the number of vertices in the graph must be equal
 * to the number of processors
 *
 * \param graph application graph structure to partition
 * \param grid grid information structure
 */
double InitialPartExec(GraphType* graph,GridType* grid){
  int i,j,v,w;
  int prov, prow;
  int istart, iend;
  int t_weight;		/**< total vertex weight */
  double p_weight;	/**< total processor power weight */
  int* glayers;
  int* vtxorder;
  int* gainorder;

  InitCtrl ICtrl; 	/**< Initial partitioning Control */
  GridType*  subg;
  Mem mem;
  MEM_Init(&mem);
  ICtrl.glayers = NULL;
  ICtrl.vtxorder = NULL;
  ICtrl.gainorder = NULL;

  int       npro   = grid->npro;
  int    nnedges   = grid->nedges;
  int   ncluster   = grid->ncluster;
  double** ratio   = grid->ratio;
  double   *exec   = grid->exec;
  double* prospeed = grid->prospeed;
  int*      rpro   = grid->rpro;
  int *cnpro = grid->cnpro;
  int*    adjidx   = graph->adjidx;
  int*    adjncy   = graph->adjncy;
  int*    adjwgt   = graph->adjwgt;
  int*      vwgt   = graph->vwgt;
  int*       pro   = graph->pro;
  int      nvtxs   = graph->nvtxs;

  int start_v;
  int sp, sv, value;

  double *iexec;
  iexec = U_malloc(npro, double, "Initial Partitioning Error Code 1\n");

  MEM_Need(&mem, 3*nvtxs, int);

  MEM_Allocate(&mem);

  glayers = ICtrl.glayers = MEM_Request(&mem, nvtxs, int);
  vtxorder = ICtrl.vtxorder = MEM_Request(&mem, nvtxs, int);
  gainorder = ICtrl.gainorder = MEM_Request(&mem, nvtxs, int);

  subg = LIBN_Init();
  subg->app = LIBG_Init();
  subg->app->nvtxs = ncluster;
  subg->app->nedges = ncluster*(ncluster-1)*2;
  LIBG_Allocate(subg->app);
  initGrid(subg, npro, nnedges);
  buildcGrid(subg, grid);

  /* First Layer Partition to different cluster */
  t_weight = 0;
  for(v = 0; v < nvtxs; v++){
    t_weight += vwgt[v];
  }

  p_weight = 0;
  for(i = 0; i < npro; i++){
    p_weight += 1/prospeed[i];
  }

  for(v = 0; v < nvtxs; v++){
    pro[v] = 0;
  }

  if(ncluster > 1){
    start_v = PseudoPeripherialVertex(0, 0, &ICtrl, graph); 
    PartExec_Domain(graph, grid, 0, ncluster-1, start_v, 0, t_weight, p_weight, &ICtrl);
  }
  for(i = 0; i < nvtxs; i++){
    if(pro[i] >= ncluster || pro[i] < 0){
      printf("vertex %d is assigned to a wrong processor %d\n", i, pro[i]);
    }
  }

  if(ncluster > 1){
    double* subexec = subg->exec;
    double* subprospeed = subg->prospeed;
    double** subratio = subg->ratio;
    int* subrpro = subg->rpro;
    int** sublatency = subg->latency;
    for(i = 0; i < npro; i++){
      subexec[i] = 0;
    }
    for(v = 0; v < nvtxs; v++){
      prov = pro[v];
      subexec[prov] += vwgt[v]*subprospeed[prov];
      istart = adjidx[v];
      iend = adjidx[v+1];
      for(i = istart; i < iend; i++){
        w = adjncy[i];
        prow = pro[w];
        if(prow != prov){
          subexec[prov] += adjwgt[i] * subratio[prov][prow];
        }
      }
    }
    for(i = 0; i < ncluster; i++){
      subrpro[i] = i;
    }
    IniProcessorCon_new(subg,graph);
    int*    sub_app_adjidx = subg->app->adjidx;
    int*    sub_app_adjncy = subg->app->adjncy;
    int*       sub_app_pro = subg->app->pro;
    for(v = 0; v < ncluster; v++){
      istart = sub_app_adjidx[v];
      iend = sub_app_adjidx[v+1];
      prov = sub_app_pro[v];
      for(i = istart; i < iend; i++){
        w = sub_app_adjncy[i];
        prow = sub_app_pro[w];
        subexec[prov] += sublatency[prov][prow];
      }
    }
    Init_Connect(graph, subg);
    InitRefine(subg, graph);
  }

  /* Second Layer Partition, partition each cluster */
  for(i = 0; i < nvtxs; i++){
    pro[i] = -pro[i];
  }
  sp = 0; /* used to compute how many processors in the cluster i*/
  for(i = 0; i < ncluster; i++){
    t_weight = 0;
    sv = -1; /* sv is the vertex belongs to cluster i */
    value = sp;
    sp += cnpro[i];
    for(j = 0; j < nvtxs; j++){
      if(pro[j] == -i){
        pro[j] = value;
        t_weight += vwgt[j];
        sv = j;
      }
    }
    start_v = PseudoPeripherialVertex(sv, value, &ICtrl, graph);
    PartExec(graph, sp-cnpro[i], sp-1, start_v, value, t_weight, &ICtrl);
  }
  
  /* Assign rpro array */
  for(i = 0; i < npro; i++){
    rpro[i] = i;
  }

  /* compute exec time */
  for(i = 0; i < npro; i++){
    exec[i] = 0;
    iexec[i] = 0;
  }
  for(v = 0; v < nvtxs; v++){
    prov = pro[v];
    exec[prov] += vwgt[v]*prospeed[prov];
    iexec[prov] += vwgt[v]*prospeed[prov];
    istart = adjidx[v];
    iend = adjidx[v+1];
    for(i = istart; i < iend; i++){
      w = adjncy[i];
      prow = pro[w];
      if(prow != prov){
        exec[prov] += adjwgt[i] * ratio[prov][prow];
        iexec[prov] += adjwgt[i] * ratio[prov][prow];
      }
    }
  }

  /****************************************/
  /* Refinement                           */
  /****************************************/
  /* copy current graph it is the same as the processor connectivity for the application */
  IniProcessorCon_new(grid,graph);

  int*    app_adjidx = grid->app->adjidx;
  int*    app_adjncy = grid->app->adjncy;
  int*       app_pro = grid->app->pro;
  int**      latency = grid->latency;
  int**     latency1 = grid->latency1;
  for(v = 0; v < npro; v++){
    istart = app_adjidx[v];
    iend = app_adjidx[v+1];
    prov = app_pro[v];
    for(i = istart; i < iend; i++){
      w = app_adjncy[i];
      prow = app_pro[w];
      iexec[prov] += latency1[prov][prow];
      exec[prov] += latency[prov][prow];
    }
  }
  double maxExec = 0;
  double maxExec1 = 0;
  for(i = 0; i < npro; i++){
    if(maxExec < exec[i]){
      maxExec = exec[i];
    }
    if(maxExec1 < iexec[i]){
      maxExec1 = iexec[i];
    }
  }

  Init_Connect(graph, grid);
  InitRefine(grid, graph);

  for(i = 0; i < npro; i++){
    exec[i] = 0;
    iexec[i] = 0;
  }
  for(v = 0; v < nvtxs; v++){
    prov = pro[v];
    exec[prov] += vwgt[v]*prospeed[prov];
    iexec[prov] += vwgt[v]*prospeed[prov];
    istart = adjidx[v];
    iend = adjidx[v+1];
    for(i = istart; i < iend; i++){
      w = adjncy[i];
      prow = pro[w];
      if(prow != prov){
        exec[prov] += adjwgt[i] * ratio[prov][prow];
        iexec[prov] += adjwgt[i] * ratio[prov][prow];
      }
    }
  }

  for(v = 0; v < npro; v++){
    istart = app_adjidx[v];
    iend = app_adjidx[v+1];
    prov = app_pro[v];
    for(i = istart; i < iend; i++){
      w = app_adjncy[i];
      prow = app_pro[w];
      exec[prov] += latency[prov][prow];
      iexec[prov] += latency1[prov][prow];
    }
  }

  maxExec = 0;
  maxExec1 = 0;
  for(i = 0; i < npro; i++){
    if(maxExec < iexec[i]){
      maxExec = iexec[i];
    }
    if(maxExec1 < exec[i]){
      maxExec1 = exec[i];
    }
  }

  MEM_Free(&mem);
  U_free("Initial Partitioning Error Code 2\n", &iexec, LEND);
  return maxExec;
}

/**
 * \Cluster level partition for each cluster 
 *
 * Partition the given graph to each cluster by its cluster compuational power
 *
 * \param graph application graph structure to partition
 * \param grid grid information structure
 * \param start_c the start cluster number
 * \param end_c the end cluster number
 * \param start_v the start vertex of the graph 
 * \param value the id presenting the region on the graph
 * \param t_weight total weight of the selected region
 * \param p_weight total weight of the processors
 * \param ctrl Initial partition control information
 */
void PartExec_Domain(GraphType* graph, GridType* grid, int start_c, int end_c, int start_v, int value, int t_weight, double p_weight, InitCtrl* ctrl){
  int i, istart, iend;
  int j;
  int v, w, prow;
  int vstart, vend;
  int nvtxs = graph->nvtxs;
  int *gadjidx = graph->adjidx;
  int *gadjncy = graph->adjncy;
  int *gvwgt = graph->vwgt;
  int *gpro = graph->pro;
  int set_value;
  int gain, changed;
  int *cnpro = grid->cnpro;
  double *cprospeed = grid->cprospeed;
  int *vtxorder = ctrl->vtxorder;
  int *gainorder = ctrl->gainorder;
  int start2, end1, start_v1=0;
  int t_limit, v_weights;
  double set_power;

  if(value != start_c){
    set_value = start_c;
  }else{
    set_value = end_c;
  }

  for(i = 0; i < nvtxs; i++){
    vtxorder[i] = -1;
    gainorder[i] = -INFINITY;
  }
  vstart = 0;
  vend = 0;

  /**< insert the gain of adjacent vertices of start_V */
  v = start_v;
  gpro[v] = set_value;
  istart = gadjidx[v];
  iend = gadjidx[v+1];
  for(i = istart; i < iend; i++){
    w = gadjncy[i];
    prow = gpro[w];
    if(prow == value){
      gain = GainCompute(w, set_value, graph);
      InsertGainOrder(gain, w, vstart, vend, vtxorder, gainorder);
      vend++;
    }
  }

  if(((end_c - start_c +1) % 2) == 0){ /* there are even number of clusters */
    set_power = 0;
    if(set_value == start_c){ /* start from the first cluster */
      for(i = start_c; i < start_c + (end_c - start_c +1)/2; i++){
        set_power += (double)cnpro[i] / cprospeed[i];
      }
    }else{ /*start from the last cluster */
      for(i = start_c + (end_c - start_c +1)/2; i < end_c+1; i++){
        set_power += (double)cnpro[i] / cprospeed[i];
      }
    }
  }else{ /* there are odd number of clusters */
    set_power = (double)cnpro[set_value] / cprospeed[set_value];
  }
  t_limit = (int) t_weight * set_power / p_weight;

  gpro[v] = set_value;
  v_weights = gvwgt[v];
  while(v_weights < t_limit){
    if(vstart < vend){
      v = vtxorder[vstart];
      vstart++;
    }else{
      for(i = 0; i < nvtxs; i++){
        if(gpro[i] == value){
          v = PseudoPeripherialVertex(i, value, ctrl, graph);
          break;
        }
      }
    }
    if(gpro[v] == value){
      gpro[v] = set_value;
      v_weights += gvwgt[v];
      istart = gadjidx[v];
      iend = gadjidx[v+1];
      for(i = istart; i < iend; i++){
        w = gadjncy[i];
        prow = gpro[w];
        if(prow == value){
          gain = GainCompute(w, set_value, graph);
          changed = 0;
          for(j = vstart; j < vend; j++){
            if(vtxorder[j] == w){
              UpdateGainOrder(gain, j, vstart, vend, vtxorder, gainorder);
              changed = 1;
              break;
            }
          }
          if(changed == 0){
            InsertGainOrder(gain, w, vstart, vend, vtxorder, gainorder);
            vend++;
          }
        }
      }
    }
  }
  
  if((end_c - start_c) > 1){
    for(v = 0; v < nvtxs; v++){
      if(gpro[v] == value){
        start_v1 = PseudoPeripherialVertex(v, value, ctrl, graph);
        break;
      }
    }
    if(((end_c - start_c +1) % 2) == 0){
      start2 = start_c + (end_c - start_c +1) / 2;
      end1 = start2 -1;
      if(value == start_c){
        PartExec_Domain(graph, grid, start_c, end1, start_v1, value, t_weight - v_weights, p_weight - set_power, ctrl);
        PartExec_Domain(graph, grid, start2, end_c, start_v, end_c, v_weights, set_power, ctrl);
      }else{
        PartExec_Domain(graph, grid, start_c, end1, start_v, start_c, v_weights, set_power, ctrl);
        PartExec_Domain(graph, grid, start2, end_c, start_v1, value, t_weight - v_weights, p_weight - set_power, ctrl);
      }
    }else{
      if(value == start_c){
        PartExec_Domain(graph, grid, start_c, end_c-1, start_v1, value, t_weight - v_weights, p_weight - set_power, ctrl);
      }else{
        PartExec_Domain(graph, grid, start_c+1, end_c, start_v1, value, t_weight - v_weights, p_weight - set_power, ctrl);
      }
    }
  }
}

/** 
 * \Inside Cluster level partition for each cluster 
 *
 * Partition the given region on each cluster
 *
 * \param graph application graph structure to partition
 * \param grid grid information structure
 * \param start_c the start cluster number
 * \param end_c the end cluster number
 * \param start_v the start vertex of the graph 
 * \param value the id presenting the region on the graph
 * \param t_weight total weight of the selected region
 * \param p_weight total weight of the processors
 * \param ctrl Initial partition control information
 */
void PartExec(GraphType* graph, int start_p, int end_p, int start_v, int value, int t_weight, InitCtrl* ctrl){
  int i, istart, iend;
  int j;
  int v, w, prow;
  int vstart, vend;
  int nvtxs = graph->nvtxs;
  int *gadjidx = graph->adjidx;
  int *gadjncy = graph->adjncy;
  int *gvwgt = graph->vwgt;
  int *gpro = graph->pro;
  int set_value;
  int gain, changed;
  int *vtxorder = ctrl->vtxorder;
  int *gainorder = ctrl->gainorder;
  int start2, end1, start_v1=0;
  int t_limit, v_weights;
  int nvtxs_value = 0;
  for(v = 0; v < nvtxs; v++){
    if(gpro[v] == value){
      nvtxs_value++;
    }
  }

  if(value != start_p){
    set_value = start_p;
  }else{
    set_value = end_p;
  }

  for(i = 0; i < nvtxs; i++){
    vtxorder[i] = -1;
    gainorder[i] = -INFINITY;
  }
  vstart = 0;
  vend = 0;

  /**< insert the gain of adjacent vertices of start_V */
  v = start_v;
  gpro[v] = set_value;
  nvtxs_value--;
  istart = gadjidx[v];
  iend = gadjidx[v+1];
  for(i = istart; i < iend; i++){
    w = gadjncy[i];
    prow = gpro[w];
    if(prow == value){
      gain = GainCompute(w, set_value, graph);
      InsertGainOrder(gain, w, vstart, vend, vtxorder, gainorder);
      vend++;
    }
  }

  if(((end_p - start_p +1) % 2) == 0){ /* there are even number of clusters */
    t_limit = (int) t_weight / 2;
  }else{ /* there are odd number of clusters */
    t_limit = (int) t_weight / (end_p - start_p +1);
  }

  v_weights = gvwgt[v];
  while(v_weights < t_limit && nvtxs_value > (int)((end_p - start_p)/2+1)){
    if(vstart < vend){
      v = vtxorder[vstart];
      vstart++;
    }else{
      for(i = 0; i < nvtxs; i++){
        if(gpro[i] == value){
          v = PseudoPeripherialVertex(i, value, ctrl, graph);
          break;
        }
      }
    }
    if(gpro[v] == value){
      gpro[v] = set_value;
      nvtxs_value--;
      v_weights += gvwgt[v];
      istart = gadjidx[v];
      iend = gadjidx[v+1];
      for(i = istart; i < iend; i++){
        w = gadjncy[i];
        prow = gpro[w];
        if(prow == value){
          gain = GainCompute(w, set_value, graph);
          changed = 0;
          for(j = vstart; j < vend; j++){
            if(vtxorder[j] == w){
              UpdateGainOrder(gain, j, vstart, vend, vtxorder, gainorder);
              changed = 1;
              break;
            }
          }
          if(changed == 0){
            InsertGainOrder(gain, w, vstart, vend, vtxorder, gainorder);
            vend++;
          }
        }
      }
    }
  }
  
  if((end_p - start_p) > 1){
    for(v = 0; v < nvtxs; v++){
      if(gpro[v] == value){
        start_v1 = PseudoPeripherialVertex(v, value, ctrl, graph);
        break;
      }
    }
    if(((end_p - start_p +1) % 2) == 0){
      start2 = start_p + (end_p - start_p +1) / 2;
      end1 = start2 -1;
      if(value == start_p){
        PartExec(graph, start_p, end1, start_v1, value, t_weight - v_weights, ctrl);
        PartExec(graph, start2, end_p, start_v, end_p, v_weights, ctrl);
      }else{
        PartExec(graph, start_p, end1, start_v, start_p, v_weights, ctrl);
        PartExec(graph, start2, end_p, start_v1, value, t_weight - v_weights, ctrl);
      }
    }else{
      if(value == start_p){
        PartExec(graph, start_p, end_p-1, start_v1, value, t_weight - v_weights, ctrl);
      }else{
        PartExec(graph, start_p+1, end_p, start_v1, value, t_weight - v_weights, ctrl);
      }
    }
  }
}

/** 
 * \Update Gain order
 *
 * when a vertex is moved into the selected region, the gain of the related vertices are changed
 * the function updates these changes
 *
 * \param gain gain value
 * \param idx the vertex 
 * \param vstart start id
 * \param vend end id
 * \param vtxorder vertex order 
 * \param gainorder gain order
 */
void UpdateGainOrder(int gain, int idx, int vstart, int vend, int* vtxorder, int* gainorder){
  int i;
  int gain1, vtx1;
  gainorder[idx] = gain;
  i = idx;
  while(i > vstart && gainorder[i] > gainorder[i-1]){
    gain1 = gainorder[i];
    vtx1 = vtxorder[i];
    gainorder[i] = gainorder[i-1];
    vtxorder[i] = vtxorder[i-1];
    gainorder[i-1] = gain1;
    vtxorder[i-1] = vtx1;
    i--;
  }
  i = idx;
  while(i+1 < vend && gainorder[i] < gainorder[i+1]){
    gain1 = gainorder[i];
    vtx1 = vtxorder[i];
    gainorder[i] = gainorder[i+1];
    vtxorder[i] = vtxorder[i+1];
    gainorder[i+1] = gain1;
    vtxorder[i+1] = vtx1;
    i++;
  }
}

/** 
 * \Insert a gain order
 *
 * when a vertex is moved into a region, there might be a new vertex gain
 * this function insert this new vertex gain
 *
 * \param gain gain value
 * \param vidx vertex id
 * \param vstart start id
 * \param vend end id
 * \param vtxorder vertex order list
 * \param gainorder gain order list
 */
void InsertGainOrder(int gain, int vidx, int vstart, int vend, int* vtxorder, int* gainorder){
  int idx;
  int gain1, vidx1;
  gainorder[vend] = gain;
  vtxorder[vend] = vidx;
  idx = vend;
  while(idx > vstart){
    if(gainorder[idx] > gainorder[idx-1]){
      gain1 = gainorder[idx];
      vidx1 = vtxorder[idx];
      gainorder[idx] = gainorder[idx-1];
      vtxorder[idx] = vtxorder[idx-1];
      gainorder[idx-1] = gain1;
      vtxorder[idx-1] = vidx1;
      idx--;
    }else{
      break;
    }
  }
}

/** 
 * \Find a pseudo Pherial Vertex 
 *
 * \param start_v start vertex of the search
 * \param value region id
 * \param ctrl Initial Partition Control information
 * \param graph the given graph
 */
int PseudoPeripherialVertex(int start_v, int value, InitCtrl* ctrl, GraphType* graph){
  int *glayers = ctrl->glayers;
  int *gpro = graph->pro;
  int *gadjidx = graph->adjidx;
  int *gadjncy = graph->adjncy;
  int nvtxs = graph->nvtxs;
  int i, istart, iend;
  int v, w;
  int morelayer, clayer;
  int return_v = start_v;

  for(v = 0; v < nvtxs; v++){
    glayers[v] = -1;
  }
  glayers[start_v] = 0;
  clayer = 0;
  morelayer = 1;
  while(morelayer == 1){
    morelayer = 0;
    for(v = 0; v < nvtxs; v++){
      if(glayers[v] == clayer){
        istart = gadjidx[v];
        iend = gadjidx[v+1];
        for(i = istart; i < iend; i++){
          w = gadjncy[i];
          if(gpro[w] == value && glayers[w] == -1){
            glayers[w] = clayer+1;
            morelayer = 1;
            return_v = w;
          }
        }
      }
    }
    clayer++;
  }
  return return_v;
}

/** 
 * \Compute the Gain
 *
 * \param v the vertex id needs to compute the gain
 * \param value the id presenting the region on the graph
 * \param graph application graph structure to partition
 */
int GainCompute(int v, int value, GraphType* graph){
  int i, istart, iend;
  int w;
  int gain;
  int *gadjidx = graph->adjidx;
  int *gadjwgt = graph->adjwgt;
  int *gadjncy = graph->adjncy;
  int *gpro = graph->pro;
  
  gain = 0;
  istart = gadjidx[v];
  iend = gadjidx[v+1];
  for(i = istart; i < iend; i++){
    w = gadjncy[i];
    if(gpro[w] == value){
      gain += gadjwgt[i];
    }else{
      gain -= gadjwgt[i];
    }
  }
  return gain;
}

/** 
 * \Local refinement, the same as the function in the uncoarsening phase
 *
 * \param n application graph structure to partition
 * \param g grid information structure
 */
void InitRefine(GridType* n, GraphType* g){
  int i, from,to,nb_moves,moves,totnb_moves,loops,nbnds;
  Mem       mem;
  Heap      pHeap,refHeap;
  RefCtrl   RCtrl;

  int     nnpro = n->npro;
  double* nexec = n->exec;

  /* prepare */
  MEM_Init(&mem);

  RCtrl.level       = -1;
  RCtrl.guide       = NULL;
  RCtrl.heap        = &refHeap;
  RCtrl.pheap       = &pHeap;
  RCtrl.g           = g;
  RCtrl.n           = n;
  RCtrl.t           = 0;

  /* alloc */
  RCtrl.procs       = U_malloc(nnpro,int,"execRefine error code 1.1\n");
  RCtrl.proref      = U_malloc(nnpro,vRefCont,"execRefine error code 1.2\n");
#if MOVE_ONCE
  RCtrl.mvtxs       = U_mallocset(RCtrl.g->nvtxs,int,0,"execRefine error code 1.3\n");
#endif
  RCtrl.bnds        = U_malloc(RCtrl.g->nvtxs,vRef,"execRefine error code 1.4\n");
  RCtrl.procon      = U_malloc(nnpro,int,"execRefine error code 1.5\n");

  nbnds             = VREF_Build(&RCtrl);

  MEM_Need(&mem,nnpro,double);
  MEM_Need(&mem,nnpro,int);
  MEM_Need(&mem,nnpro,int);
  MEM_Need(&mem,MAX_KLMOVES,vRef);
  MEM_Need(&mem,8*nbnds/2,vRef*);
  MEM_Need(&mem,nnpro,double);
  HeapInit(&mem,&pHeap,nnpro);
  HeapInit(&mem,&refHeap,8*nbnds/2);

  MEM_Allocate(&mem);

  RCtrl.flow        = MEM_Request(&mem,nnpro,double);
  RCtrl.procs2      = MEM_Request(&mem,nnpro,int);
  RCtrl.procon2     = MEM_Request(&mem,nnpro,int);
  RCtrl.undo        = MEM_Request(&mem,MAX_KLMOVES,vRef);
  RCtrl.nb_undo     = 0;
  RCtrl.exec        = MEM_Request(&mem,nnpro,double);
  RCtrl.nrvtxs      = 0;
  RCtrl.rvtxs       = MEM_Request(&mem,8*nbnds/2,vRef*);

  for(i=0;i<nnpro;i++)
    RCtrl.t += nexec[i];

  RCtrl.tmean = RCtrl.t/nnpro;
  RCtrl.tmax = 0;
  RCtrl.tmin = INFINITY;

  for(i=0;i<nnpro;i++){
    if(nexec[i] > RCtrl.tmax)
      RCtrl.tmax = nexec[i];
    if(nexec[i] < RCtrl.tmin)
      RCtrl.tmin = nexec[i];
  }

  RCtrl.im          = nnpro*RCtrl.tmax/RCtrl.t;
  RCtrl.imthreshold = IMTHRESHOLD(MAX_IM,levels,&RCtrl);

  RCtrl.gain_fn   = &OneGain;

  totnb_moves = 0;
  nb_moves    = 1;
  loops = 0;

  while(nb_moves != 0){
    nb_moves = 0;

    /* from overloaded */
    RCtrl.filter_fn = &OneFromFilter;
    HeapClean(&pHeap);
    for(i=0;i<nnpro;i++)
      HeapInsert(&pHeap,nexec[i]);

    while((from = HeapGetHead(&pHeap)) != -1){
      if(!FILTEREXECFROM(&RCtrl,nexec[from]))
        break;
      HeapifyFrom(&RCtrl,from);
      if(HeapGetHead(RCtrl.heap) != -1){
        if((moves = IRefine(&RCtrl)) != 0){
          nb_moves += moves;
        }
      }else{
        HeapRemoveHead(&pHeap);
      }
    }

    /* to underloaded */
    RCtrl.filter_fn = &OneToFilter;

    HeapClean(&pHeap);

    for(i=0;i<nnpro;i++)
      HeapInsert(&pHeap,-nexec[i]);

    while((to = HeapGetHead(&pHeap)) != -1){
      if(!FILTEREXECTO(&RCtrl,nexec[to]))
        break;

      HeapifyTo(&RCtrl,to);
      if(HeapGetHead(RCtrl.heap) != -1){
        if((moves = IRefine(&RCtrl)) != 0){
          nb_moves += moves;
        }
      }else{
        HeapRemoveHead(&pHeap);
      }
    }
    totnb_moves += nb_moves;
  }
  /******************************************************************************************/
#if DIFFUSE
  for(i=0;i<nnpro;i++)
    RCtrl.flow[i] = nexec[i] - RCtrl.tmean;

  totnb_moves = 0;

  if(RCtrl.im>RCtrl.imthreshold){
    RCtrl.gain_fn   = &TwoGain;

    nb_moves    = 1;
    loops = 0;

    while(nb_moves != 0){
      nb_moves = 0;

      /* from overloaded */
      RCtrl.filter_fn = &TwoFromFilter;

      HeapClean(&pHeap);

      for(i=0;i<nnpro;i++)
        HeapInsert(&pHeap,nexec[i]);

      while((from = HeapGetHead(&pHeap)) != -1){
        if(!FILTEREXECFROM(&RCtrl,nexec[from]))
          break;
        HeapifyFrom(&RCtrl,from);
        if(HeapGetHead(RCtrl.heap) != -1){
          if((moves = IRefine(&RCtrl)) != 0){
            nb_moves += moves;
          }
        }else
          HeapRemoveHead(&pHeap);
      }

      /* to underloaded */
      RCtrl.filter_fn = &TwoToFilter;

      HeapClean(&pHeap);

      for(i=0;i<nnpro;i++)
        HeapInsert(&pHeap,-nexec[i]);

      while((to = HeapGetHead(&pHeap)) != -1){
        if(!FILTEREXECTO(&RCtrl,nexec[to]))
          break;

        HeapifyTo(&RCtrl,to);
        if(HeapGetHead(RCtrl.heap) != -1){
          if((moves = IRefine(&RCtrl)) != 0){
            nb_moves += moves;
          }
        }else
          HeapRemoveHead(&pHeap);

        if(RCtrl.im<RCtrl.imthreshold)
          break;
      }

      totnb_moves += nb_moves;
    }
  }

#endif
  debugExecNoPrint(&RCtrl);
  VREF_Free(&RCtrl);
  U_free("execRefine error code 2.1\n",&(RCtrl.procs),&(RCtrl.proref),&(RCtrl.bnds),LEND);
#if MOVE_ONCE
  U_free("execRefine error code 2.2\n",&(RCtrl.mvtxs),LEND);
#endif

  MEM_Free(&mem);
}

/** 
 * \Refinement process, the same as the function in the uncoarsening phase 
 *
 * \param ctrl Initial Refinement control information
 */
int IRefine(RefCtrl* ctrl){
  int v,idx,from,to,nb_moves = 0,dOflag;
  double deltaF,deltaT,deltaO;

  FILTER_FN* filter_fn = ctrl->filter_fn;

  Heap*  heap          = ctrl->heap;
  vRef** rvtxs         = ctrl->rvtxs;
  vRef*  rvtx;

  while((v = HeapRemoveHead(heap)) != -1){
    rvtx   = rvtxs[v];
    idx    = rvtx->idx;
    to     = rvtx->to;
    from   = rvtx->from;
    deltaF = rvtx->deltaF;
    deltaT = rvtx->deltaT;
    deltaO = rvtx->deltaO;
    dOflag = rvtx->dOflag;

    /********************************/
    /* reassignment passing filters */
    /********************************/
    if((* filter_fn)(ctrl,1,idx,from,to,deltaF,deltaT,deltaO,dOflag)){
      /* update rinfo was deleted after that */
      MOVE_NOUNDO(ctrl,rvtx);
      nb_moves++;
      continue;
    }

    /*************************************/
    /* not an interesting move remove it */
    /*************************************/
    HeapRemove(heap,v);
    AL_REM(rvtxs + v,rvtxs,ctrl->nrvtxs);
  }
  return nb_moves;
}

/** 
 * \build a cluster level grid
 *
 * build a grid treating each cluster in the given grid as a processor
 *
 * \param subg the cluster level grid
 * \param g grid information structure
 */
void buildcGrid(GridType* subg, GridType* g){
  int* subadjidx = subg->adjidx;
  int* subadjncy = subg->adjncy;
  int* subadjwgt = subg->adjwgt;
  int* subadjlcy = subg->adjlcy;
  int** sublatency = subg->latency;
  double* subprospeed = subg->prospeed;
  double** subratio = subg->ratio;

  int* cnpro = g->cnpro;
  double* cprospeed = g->cprospeed;
  int* cadjidx = g->cadjidx;
  int* cadjncy = g->cadjncy;
  int* cadjwgt = g->cadjwgt;
  int* cadjlatency = g->cadjlatency;
  int* cadjfrom = g->cadjfrom;
  int* cadjto = g->cadjto;
  double** ratio = g->ratio;

  int ncluster = g->ncluster;
  int nnedges = g->nedges;

  int i, j, istart, iend;
/* build a cluster based grid */
  subg->npro = ncluster;
  double max_pro = 0;
  for(i = 0; i < ncluster; i++){
    if((cnpro[i] / cprospeed[i]) > max_pro){
      max_pro = cnpro[i] / cprospeed[i];
    }
  }
  for(i = 0; i < ncluster; i++){
    subprospeed[i] = max_pro / (cnpro[i] / cprospeed[i]);
  }
  for(i = 0; i < ncluster+1; i++){
    subadjidx[i] = cadjidx[i];
  }
  for(i = 0; i < nnedges; i++){
    subadjncy[i] = cadjncy[i];
    subadjwgt[i] = cadjwgt[i];
    subadjlcy[i] = cadjlatency[i];
  }
  for(i = 0; i < ncluster; i++){
    for(j = 0; j < ncluster; j++){
      sublatency[i][j] = INFINITY;
    }
    sublatency[i][i] = 0;
  }
  for(j = 0; j < ncluster; j++){
    istart = subadjidx[j];
    iend = subadjidx[j+1];
    for(i = istart; i < iend; i++){
      sublatency[j][subadjncy[i]] = subadjlcy[i];
      subratio[j][subadjncy[i]] = ratio[cadjfrom[j]][cadjto[j]];
    }
  }
}
